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Abstract. A new computational analysis tool, downscaling test, is introduced and applied for studying the 
convergence rates of truncation and discretization errors of finite-volume discretization schemes on general irregular 
(e.g., unstructured) grids. The study shows that the design-order convergence of discretization errors can be achieved 
even when truncation errors exhibit a lower-order convergence or, in some cases, do not converge at all. The down- 
scaling test is a general, efficient, accurate, and practical tool, enabling straightforward extension of verification 
and validation to general unstructured grid formulations. It also allows separate analysis of the interior, boundaries, 
and singularities that could be useful even in structured-grid settings. There are several new findings arising from 
the use of the downscaling test analysis. It is shown that the discretization accuracy of a common node-centered 
finite- volume scheme, known to be second-order accurate for inviscid equations on triangular grids, degenerates to 
first order for mixed grids. Alternative node-centered schemes are presented and demonstrated to provide second 
and third order accuracies on general mixed grids. The local accuracy deterioration at intersections of tangency 
and inflow/outflow boundaries is demonstrated using the DS tests tailored to examining the local behavior of the 
boundary conditions. The discretization-error order reduction within inviscid stagnation regions is demonstrated. 
The accuracy deterioration is local, affecting mainly the velocity components, but applies to any order scheme. 
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Regular grids 

• Structured mapped grids 


Irregular grids 

• Structured nonmapped grids 

• Unstructured grids 


Fig. 1.1. Grid classification. 


1. Introduction. There is an increasing reliance on computational simulations in air- 
craft design practices, supplementing traditional analytic and experimental approaches. Veri- 
fication and validation methodologies [20] are being developed to ensure the correct applica- 
bility of these approaches in practical applications on various grids. The following classifica- 
tion is used to characterize interior grids: grids with periodic (repeating) topological structure 
are called structured ; grids derived by smooth mapping from structured grids with periodic 
metrics are called mapped grids. They are referred in this paper as regular grids (examples 
include, but are not limited to, grids derived from Cartesian ones: triangular grids obtained by 
the same diagonal splitting, stretched grids, smooth curvilinear grids, etc.). Grids with vary- 
ing local topology are called unstructured , e.g., grids with the number of neighboring nodes 
changing from node to node; unstructured grids and nonmapped structured grids are consid- 
ered irregular. Figure 1.1 illustrates this grid classification. Verification methodologies for 
regular grids, e.g., [17], are relatively well-developed in comparison to irregular grids, espe- 
cially unstructured grids containing mixed elements or grids derived through agglomeration 
techniques. The summary of the latest of three Drag Prediction Workshops [14] illustrates the 
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problems associated with assessing errors in practical complex-geometry/complex-physics 
applications. 

The two-dimensional (2D) primal meshes considered in this paper are composed of tri- 
angular and quadrilateral cells. The finite-volume discretization (FVD) schemes are derived 
from the integral form of a steady conservation law 

/(F-fi) dT = JJ{f-S) dfi, (1.1) 

r n 

where / is a forcing function independent of the solution, S is a solution-dependent source 
function, fl is a control volume with boundary T, n is the outward unit normal vector, and 
F is the flux vector that may include viscous and/or inviscid contributions. The general 
FVD approach requires partitioning the domain into a set of non-overlapping control volumes 
and implementing numerically equation (1.1) over each control volume. The focus of this 
paper is node-centered FVD schemes, where solution values are stored at the mesh nodes; 
the proposed analysis techniques have already been applied in [25] to cell-centered FVD 
schemes. 



FIG. 1.2. Median-dual partition for node-centered finite-volume discretizations. Po — P 4 denote grid nodes. 
The control volume (dual cell ) around Pq is shaded. 


For 2D node-centered FVD schemes, a median-dual partition can be constructed by con- 
necting the mass centers of the primal-mesh cells with the midpoints of the surrounding edges 
(Figure 1.2). These non-overlapping control volumes cover the entire computational domain 
and compose a mesh that is dual to the primal mesh. 

The discrete solution is represented as a piecewise polynomial function; the polyno- 
mials are defined either within primal or dual cells. With discrete solution defined at the 
grid nodes, the fluxes at the dual-cell boundaries are reconstructed using these local polyno- 
mial representations. The discretization is applied at a sequence of refined grids satisfying 
the consistent refi nement property. The property requires the maximum distance across pri- 
mal and dual cells to decrease consistently with increase of the total number of grid points, 
N. In particular, the maximum distance should tend to zero as N^ 1 ^ 2 in 2D computations. 
For three-dimensional (3D) unstructured grids, the consistent refinement property is studied 
in [25]. 

The main accuracy measure of any FVD scheme is the discretization error , Ed, defined 
as the difference between the exact discrete solution, Q h , of the discretized equations (1.1) 
and the exact continuous solution, Q, to the differential conservation law 

VF = f — S, 
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E d = Q — Q h ; 


(1.3) 


Q is sampled at grid nodes. 

A common approach to evaluate the accuracy of a FVD scheme is to monitor the con- 
vergence of truncation errors. Truncation error, E t , measures the accuracy of the discrete 
approximation to the differential equations ( 1 .2). For finite differences, it is defined as the 
residual obtained after substituting the exact solution Q into the discretized differential equa- 
tions [15], For FVD schemes, the traditional truncation error is usually defined from the 
time-dependent standpoint [23, 24], In the steady-state limit, it is defined (e.g., in [10]) as the 
residual computed after substituting Q into the normalized discrete equations (1.1), 



JJ ( f h - S h ( q )) cm + J (F h (Q) ■ ri) dr , 


n 


r 


(1.4) 


where F h is a reconstruction of the flux F at the boundary F, |0| is the measure of the control 
volume. 


|f2| = JJ dQ, (1.5) 

n 

f h and S h are, respectively, approximations of the forcing function / and the source function 
S on Q, and the integrals are computed according to some quadrature formulas. Conver- 
gence of truncation errors has been applied as a FVD accuracy measure on both regular and 
irregular grids [5, 22], On regular grids, this approach is well justified; the choice of scaling 
factor y^y ensures that truncation errors converge as 0(h p ) on sequences of consistently re- 
fined grids, where h is a characteristic meshsize and p is the design discretization-accuracy 
order of the method. However, the truncation-error convergence is often misleading for FVD 
schemes defined on irregular (e.g., unstructured) grids. Several studies, e.g., [19, 22], noted 
that 2" ' / -order convergence of truncation errors for some commonly used node-centered FVD 
schemes can be achieved only on grids with a certain degree of geometric regularity. Other 
studies, e.g., [10, 13], showed that such degradation of truncation-error convergence does 
not always imply a lower-order convergence of discretization errors. Examples in subsequent 
sections confirm that on irregular grids, the design-order discretization-error convergence can 
be achieved even when truncation errors exhibit a lower-order convergence or, in some cases, 
do not converge at all. Note that these results do not contradict to the Lax theorem, which 
states that consistency (convergence of truncation errors) and stability are sufficient (not nec- 
essary) for convergence of discretization errors. In fact, for some formally inconsistent FVD 
schemes, it has been rigorously proven that the discretization errors converge [7]. 

Although the convergence of irregular-grid truncation errors is not identical to the discre- 
tization-error convergence, it still can be monitored to indicate if the design order of the 
discretization accuracy can be achieved. On a sequence of irregular grids satisfying the con- 
sistent refinement property, the convergence order of truncation errors is typically less than 
the discretization-error convergence order by 1 for inviscid equations and by 2 for viscous 
equations. 

The main computational tool introduced in this paper for evaluating accuracy of dis- 
cretization schemes is a downscaling (DS) test. The test is performed for a known exact 
or manufactured solution. For a manufactured solution, the forcing functions are found by 
substituting this solution into the continuous governing equations and boundary conditions. 
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The DS test consists of a series of inexpensive computational experiments that account for 
local properties of the studied scheme; it is designed to provide estimates for the convergence 
orders of the discretization and truncation errors by comparing errors obtained on differ- 
ent scales. The DS test is a very general technique that can be applied to arbitrary grids and 
geometries. It can be tailored to study the discretization accuracy in the interior, at the bound- 
ary, and/or in vicinity of singularities. Analysis methods predicting the performance of DS 
tests have also been developed. Although, the computations presented in this paper are one- 
and two-dimensional, the analysis techniques and conclusions are fully extensible and have 
already been successfully applied in the three dimensions; some relevant 3D computations 
have been reported in [25]. While rigorous proofs of estimates for convergence of FVD trun- 
cation and discretization errors in realistic computations are still out of reach, the introduced 
practical analysis methods enable accurate and efficient accuracy assessment on various ir- 
regular grids, for various types of solutions and boundary conditions. Numerous examples 
reported here and in other recent publications ([8, 25]) confirm the predicting power and wide 
applicability of the DS-analysis. 

The material in this paper is organized as follows. Section 2 describes the construction 
and application of DS tests for general FVD schemes on irregular grids. Section 3 introduces 
analytical methods for predicting convergence of discretization and truncation errors in a DS 
test. Section 4 includes one-dimensional examples providing insights into convergence of the 
discretization and truncation errors on irregular grids. Section 5 analyzes convergence for 
several sets of 2D equations and FVD schemes and corroborates the analysis with numerical 
tests performed on irregular triangular, quadrilateral, and mixed-element grids. In addition, 
the effects of flow singularities are analyzed; one-order reduction of the discretization ac- 
curacy is observed and explained for inviscid stagnation flows. The concluding Section 6 
summarizes the demonstrated computational and analytical results and discusses complexity 
of the introduced FVD schemes. 

2. Downscaling (DS) test. The purpose of the downscaling (DS) test is to predict the 

discretization and truncation-error convergence orders in computations performed on irregu- 

lar grids satisfying the consistent refinement property. To apply a DS test, one first chooses a 
representative manufactured solution defined on the given computational domain. 

The DS test requires numerical computations on a sequence of contracted domains zoom- 
ing toward a focal point within the original computational domain (Figure 2.1). The choice of 
the focal-point location can be varied to study a typical interior discretization, a boundary dis- 
cretization, or a particular singularity. Because the number of points in the DS-test domains 
is held (approximately) fixed, one can study solutions on grids with characteristic meshsizes 
much smaller than those affordable in computations with a globally-refined grid sequence. 

There are at least two possible strategies for grid generation on these contracted domains: 

(1) The first strategy is termed “scaled grid”; with this strategy, the first (coarsest) computa- 
tional domain is defined as a subdomain of the investigated global mesh containing the focal 
point; other (finer) domains and their mesh patterns are derived by scaling down this first do- 
main (e.g., repeatedly multiplying all the distances from the focal point by a given factor, say, 

1/2 or 1/4). (2) In an alternative, independent-grid approach, a new grid can be generated on 
each domain, assuming that the consistent refinement property is satisfied, i.e., the maximum 
distance across a grid cell is scaled down with the same rate as the diameter of the contracted 
domains. 

The “scaled-grid 2 * * * * * * * * 11 approach (sketch (a) in Figure 2.1) is especially useful for studying 
interior discretizations and straight boundaries. It is impractical for studies near a general 

(discretely defined) curvilinear boundary because the physical boundary shape should be pre- 

served on each scale in the DS sequence (sketch (c) in Figure 2.1). In the DS test, the studied 
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(a) Scaled-grid approach 


(b) Independent-grid ap- 
proach 


(c) Accounting for curved 
physical boundary 


FIG. 2.1. DS test: illustration of computational domains for studying the interior and general boundary. Black 
bullets mark the focal points; white squares mark the interface between the interior and the DS-test domain. 


FVD scheme is supplemented with a set of boundary conditions at the interface between 
the interior and the DS domain; overspecification from the known manufactured solution 
is a usual choice. The freedom to choose the focal point, the shape of domains, the type 
of boundary conditions, and the grid generation method greatly simplifies DS testing. In 
studying accuracy of discrete boundary conditions, i.e., placing the focal point at the physical 
boundary, the physical conditions are implemented at the studied boundary surface; overspec- 
ification can be applied at the interface between the interior and the DS domain (see sketches 
in Figure 2.1). 

The DS test evaluates local truncation and discretization-error convergence orders by 
comparing errors obtained in computations on different scales. The DS convergence of trun- 
cation errors is an accurate indication of the truncation-error convergence (in the Toe- norm) 
observed in global grid-refinement computations, providing the DS test samples all represen- 
tative regions. Convergence of grid-refinement truncation errors measured in integral norms, 
e.g., L-\ -norm, may be higher order because these norms are less sensitive to fluctuations 
occurring locally. 

On grids with an increased regularity, the DS test may overestimate convergence of the 
discretization errors. Some global error accumulation may occur on such grids; because 
of its local nature and overspecified interface boundary conditions, the DS test is incapable 
to account for the global error accumulation. Our experience shows that on truly irregular 
multidimensional grids, errors produced locally dominate globally accumulated errors, and 
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the discretization-error estimates obtained in a DS test become sharp. 

3. Analysis of discretization and truncation errors. In this section, we describe an 
analytical approach to estimating the convergence orders of the discretization and truncation 
errors on irregular grids. 

3.1. Foundations of analysis. The analysis is performed for a known exact or manu- 
factured solution. For a given FVD scheme, we define a general residual function 


R{q) 


1 

P 


/ 


,r 


(F h (q) ■ n) dr 


JJ if h - S h (q)) dfi 

n 


(3.1) 


where q is a discrete function defined at the grid nodes; j3 is a scaling factor; at each control 
volume, the flux reconstruction F h and the source approximation S h are computed for q, 
using local polynomial approximation; and the outward unit normal n, the discrete force 
function f h , and all integrals are approximated with certain accuracy by predefined methods. 
The general residual accounts for interior discretization and all boundary conditions. Note 
that by definition 


R(Q h ) = 0 (3.2) 

and 

R(Q ) = (3.3) 

Recall that Q h and Q are the exact discrete solution and a discrete (sampled-at-nodes) repre- 
sentation of the exact continuous solution, respectively. Substituting (1.3) into (3.2), 


R(Q - E d ) = 0, (3.4) 

and assuming the discretization error to be small comparing to the exact continuous solution, 
Q, (\E d \ « |Q|), we obtain 


R{Q) - J(Q)E d » 0, 

where J (Q) is the Jacobian of R(q) computed for q = Q. 


J (Q) = 


18 _ 

P dq 


j (F h (Q) -h)dT + jj S h (Q) dfl 


Lr n 

Thus, the discretization error can be estimated as 

E d ~ J~ 1 (Q)R(Q), 


(3.5) 


(3.6) 


(3.7) 


The equation (3.7) is the main vehicle for relating discretization and truncation errors 
and for estimating the discretization-error convergence order. Note that the estimate does not 
depend on the scaling factor (3, so 0 can be adjusted to simplify analysis. In particular, in the 
course of this paper, we choose 


P= |r|, 


6 


(3.8) 



where |T| = 1 ^ d = ft d_1 is the measure of the control-volume boundary, d is the 

space dimension, and ft, is a characteristic diameter of the control volumes. Note, that with 
this scaling, the truncation error computed according to (1.4) relates to R(Q) as 



Estimates of residual (and truncation-error) convergence are relatively easy to obtain for 
a given manufactured solution because the estimation does not require solution of discrete 
equations; residuals can be directly measured in computations. Thus, the main complexity of 
evaluation of the discretization-error convergence order rests with evaluation of the inverse 
Jacobian, which accounts for both interior and boundary discretizations. Example evaluations 
of the inverse Jacobian for formulations focusing on boundary conditions are provided in [9]. 

3.2. Local analysis of Jacobian . The asymptotic order of the inverse Jacobian can 
be evaluated locally. For general systems of nonlinear equations, the asymptotic order, mj, 
J~ 1 (Q) = 0(h mj ) can be predicted by analyzing an equivalent linear operator , E(Q). The 
equivalent linear operator is derived locally by a scaling analysis of the J(Q) operator (3.6); 
the analysis evaluates magnitudes of the coefficients appearing in the matrix of the discrete 
Jacobian. In the scaling analysis, all discrete spatial derivatives operators contributing to the 
Jacobian are replaced with 0(l/ft fe ) terms, where k is the differentiation order. Surface inte- 
grals, i.e., terms like f (g{Q)n x ) dT, are replaced with g(Q)0(h d ~ 1 ) and volume integrals, 
r 

i.e., terms like J j g(Q) dil, are replaced with g(Q)0(h d )\ here g(Q) is a functional of the 

n 

(manufactured) solution and/or its derivatives and n x is the .r-dircctional component of the 
outward unit normal n. We will refer to these replacements as equivalent substitutions. After 
equivalent substitutions, E(Q) can be formally inverted to estimate the order of ■J~ 1 {Q). For 
basic fluid equations, this analytical evaluation can be performed easily, as shown in subse- 
quent examples. 

Some simple practical rules can be suggested for quick evaluation of m j for non-degene- 
rate equations. A scalar (nonlinear) equation is called non-degenerate for a given manufac- 
tured solution, if the coefficients of the highest derivatives in the equation linearized around 
this solution do not all vanish; a system of differential equations is non-degenerate if the de- 
terminant of its linearization is non-degenerate. Note that change of variables does not affect 
equation degeneration; it is determined by the type of equation and the manufactured solution. 
With the chosen scaling factor (3.8), integration over dual-cell boundaries does not introduce 
an h-dependent factor to the Jacobian; control volume integration introduces a factor O(h). 

For FVD schemes discretizing non-degenerate scalar conservation laws with zero source 
function S, 


m.j = (order of solution differentiation in F). (3.10) 

For example, for inviscid fluxes, mj = 0; for viscous fluxes, m.j = 1. For system of equa- 
tions, the operator J _1 (Q) is a matrix, and the asymptotic order mj may vary for different 
matrix entries. Examples of evaluation the asymptotic order for systems of equations are 
given subsequently for non-degenerate equations in Section 5.2.1 and for degenerate equa- 
tions in Section 5.3 . 

3.3. Error accumulation. The assessment of sources for possible error accumulation 
presented in this section follows the arguments provided in [10]. The discretization error can 
be considered as a convolution of the truncation-error function with the corresponding Green 
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function. Assuming that the Green function is smooth almost everywhere, there are two pos- 
sible scenarios: (1) There is no significant truncation-error cancellation in the convolution 
integration. This scenario corresponds to global error accumulation and produces discretiza- 
tion errors converging with the same order as truncation errors. (2) There is a significant 
truncation-error cancellation in the convolution integration. In this case, locally produced 
errors dominate globally accumulated errors, and discretization errors converge with a higher 
order than truncation errors. 

From (3.1), there are two contributors to the residual of an FVD scheme: flux computa- 
tions and source/force computations. On irregular grids, the flux-related contributions do not 
cause dominant global discretization-error accumulation. The rationale is that a flux error at 
an internal face contributes equal and opposite amounts to the residuals on both sides of the 
face; thus, the generated error is (almost) canceled in the residual integration with a smooth 
convolution kernel. Some error accumulation may still occur because of variations in the con- 
volution kernel, which is smooth, but not constant. On regular grids, this accumulated error 
might become dominant because the leading contributions to the local residuals are canceled 
out. The net result is that the discretization-error convergence order is expected to be the 
same on regular and irregular grids. 

The source and force residual contributions at a control volume are independent of neigh- 
boring control volume residuals. Large smooth residual components can be introduced by the 
source/force integration method. To prevent accumulation of a dominant error, the integration 
method should provide a design-order accuracy for integration of a smooth function over any 
(){ 1 ) -large domain. Note that the local integration accuracy on individual control volumes is 
allowed to be lower order. An example is shown in Section 5.1. 

3.4. Analysis of residual convergence. While residuals can be directly measured in 
computations, for scheme design purposes, it is useful to analyze separate factors contributing 
to the residual and determine the asymptotic residual order, niR, R(Q) = 0(h mR ). The 
asymptotic order is evaluated for the chosen scaling factor (3.8). 

3.4.1. Analysis of residual: flux contributions. There are three main flux-related con- 
tributors to R(Q): (1) flux reconstruction, (2) dual-cell boundary approximation, and (3) flux 
integration; these contributions have also been identified in [18]. A p th - order residual R(Q) 
(mu = p ) can be obtained, if the following three sufficient accuracy conditions are satisfied. 
The conditions are formulated in terms of the maximal degree of a polynomial, for which the 
considered approximation provides the exact outcome. 

1. Flux reconstruction accuracy. The flux F h should be reconstructed precisely at the 
control-volume boundaries for an analytical flux F represented by polynomials of the (p — 
l) th degree. 

2. Control-volume boundary accuracy. The integral of the outward normal to the bound- 
ary should be computed precisely for each boundary segment (face) represented as a (p— l) th - 
order curve, i.e., its (piecewise) parametric formulation with respect to the curve length in- 
volves only polynomials of the (p— l) th degree. In many practical cases, the control volumes 
in the interior have piecewise linear boundaries, so this contributor plays an important role 
only near the physical boundaries of the computational domain, where the imposed geometry 
must be approximated. 

3. Flux integration accuracy. The flux integration method should be exact for fluxes 
represented as polynomials of the (p — l) th degree. Note that this condition relates to the 
accuracy of the integration scheme, and assumes no reconstruction errors. In particular, in 
ID FVD schemes, this condition is satisfied for all polynomials because each of the two 
control-volume boundaries collapse to a single point. In multiple dimensions, one should 



distinguish between the global and local integration accuracies. The global p th -order integra- 
tion accuracy is achieved when the numerical integration of polynomials of {p — l) t/l degree 
over the entire closed dual-cell boundary is exact, i.e., provides the same result as the analyt- 
ical integration. To achieve the local integration accuracy of p th order, the (p — l) th degree 
polynomials should be integrated exactly over each segment (face) of the boundary. Local 
accuracy is more strict and implies global accuracy of the same order or higher; global p th - 
order flux integration accuracy is sufficient. For many FVD schemes, local and global flux 
integration accuracies are the same; however, one commonly used FVD scheme discussed in 
Section 5.2.2 relies on an economic integration strategy that capitalizes on the global integra- 
tion accuracy order exceeding the order of the local integration accuracy. 

Each of these sufficient accuracy conditions can be evaluated separately and indepen- 
dently. Note: (1) the flux reconstruction and flux integration conditions recover the steps 
given in [4] for linear and quadratic schemes on unstructured meshes; (2) these conditions are 
conservative (not necessary) because higher-order approximation to R I Q ) can be achieved 
due to error cancellation, especially, on regular grids. 

3.4.2. Analysis of residual: source and force contributions. Generally, to provide an 
0{h p ) contribution to the residual R ( Q ), it is sufficient for the numerical integration over 
the control volume to be exact for the integrated function represented by a polynomial of 
the (p — 2) th degree. However, because of possible global error accumulation, this is not 
sufficient; the design -order accuracy should be demonstrated in integration of an arbi- 
trary smooth function over any ()( l)-large domain. To verify an integration scheme, one can 
(1) choose a continuous function and integrate it analytically over an 0(l)-large computa- 
tional domain; (2) implement the studied numerical integration scheme over this domain on 
a sequence of globally refined meshes; and (3) monitor convergence of the integration error, 
which is expected to be 0{h p ), where p is the design order. 

3.5. Error estimation and applications. The methodology for error estimation pro- 
posed in this paper suggests two ways to evaluate asymptotic convergence of truncation and 
discretization errors that would be observed in global grid-refinement computations on irreg- 
ular grids. 

The main computational analysis tool is a DS test that explicitly demonstrates error con- 
vergence which is used as predictor for the error convergence in global computations. The 
predictions of truncation-error convergence are always accurate, providing sufficient sam- 
pling of DS computational windows (i.e., if DS tests are performed in all representative lo- 
cations). Predictions of discretization-error convergence may not be accurate on grids with 
increased regularity. DS tests are especially useful for scheme verification. Slow convergence 
of discretization errors detected in a DS test is an unambiguous indication of deficiency of 
the analyzed scheme. 

Another, pure analytical, method is to use estimates for asymptotic orders of the inverse 
Jacobian and the residual. The analysis is especially useful for FVD scheme design because 
it identifies error contributions from separate parts of the scheme. For conservation laws 
with J -1 (Q) = 0(h mj ) and R ( Q ) = 0(h mR ), the truncation and discretization errors are 
evaluated from (3.9) and (3.7), respectively, as 

E t = 0(h mR ~ 1 ) and E d = 0(h mR+mj ). (3.11) 

Both the DS tests and the estimates (3.11) should be complimented with a verification of 
the source/force integration scheme for global computations. 

4. One-dimensional examples. Examples in this section show FVD schemes on one- 
dimensional (ID) irregular grids and illustrate applications of the analysis and DS tests to 
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predicting error convergence. For ID equations, sufficient conditions for control-volume 
boundary accuracy and for flux-integration accuracy (conditions 2 and 3 in Section 3.4.1) are 
automatically satisfied. For simplicity, the examples considered in this section do not have 
source terms and the forcing term is integrated analytically, so itir is fully determined by the 
flux reconstruction accuracy. 

4.1. Discretization grids. A ID discretization grid is defined as a combination of the 
primal and dual nodes. The solution values are located at the primal nodes; the fluxes are 
located at the dual nodes. For node-centered discretizations, a natural strategy is to place the 
primal mesh first and, then, use this mesh as a reference for placing dual control volumes. 

The ID node-centered discretization grids employed in this section are designed to study 
effects of grid irregularities and are described as follows. The first and the last of the N + 
1 primal nodes, X{, i = 0,1 , ... ,7V, are always located at the ends of the computational 
interval; the interior nodes can be distributed either uniformly or randomly. Either distribution 
retains the nodal ordering and ensures that the maximal distance between the neighboring 
nodes is 0(1/N). Let s*, i = 0, 1, . . . , (N + 1) denote the flux locations. The first and the 
last fluxes are also located at the ends of the interval; the location of an interior flux, S j, is 
always between the primal nodes Xi-i and x, and initially defined as either a biased or an 
unbiased average; then, the dual node S j may be randomly perturbed. 


+ I I • - 




(a) Uniform primal mesh; unperturbed unbiased dual mesh 


♦-+•1 


i • 


(b) Uniform primal mesh; perturbed unbiased dual mesh 
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(c) Random primal mesh; unperturbed unbiased dual mesh 


FIG. 4.1. Examples of one-dimensional discretization grids: black bullets denote primal mesh nodes, vertical 
tic-marks denote dual mesh nodes 


Specifically, on an interval x £ [a, 6], the primal nodes are distributed according to 

x 0 = a; Xi = a + (i + r,) *= 1, . . . , (iV - 1); x N =b; (4.1) 

where Vi is either zero (uniform primal mesh) or a random number —0.4 < r, < 0.4 (random 
primal mesh). The dual-mesh nodes are computed accordingly as 

s 0 = a; Si = Xi-i + di (xi - Xi-i) , i = s N+ i = b; (4.2) 

where di = 0.5 corresponds to an unbiased unperturbed dual mesh; di = 0.7 corresponds to a 
biased unperturbed dual mesh; and di = 0.5 + rf or di = 0.7 + rf correspond to unbiased and 
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biased perturbed dual meshes, respectively; here r\ is a random number —0.25 < rf < 0.25. 
Grid examples are shown in Figure 4. 1 . Note, that a combination of a random primal mesh 
and an unbiased, perturbed dual mesh is representative for multidimensional median-dual 
partition on irregular grids. 

The global computations in this section refer to tests performed on the interval x £ 
[0, 1] using a sequence of grids with the total number of grid nodes increasing as N = 
2 3 , 2 4 , . . . , 2 14 . The DS test is performed on a sequence of the intervals [0.5 — /, 0.5 + Z], 
where l = 2 _1 , 2~ 2 , . . . , 2 -10 is a scaling factor; the random grids on each scale are 
generated independently; the number of grid nodes at each scale is fixed, N = 16. 

4.2. Convection equation. In the first example concerned with a constant-coefficient 
convection equation, we illustrate applications of the DS test to predict the accuracy of irreg- 
ular FVD schemes. The equation 

d x u = f{x), 17(0) = U 0 (4.3) 

is satisfied with the exact solution U = sin (x) , f = cos (x). This is a constant-coefficient 
equation and, therefore, non-degenerate. FVD equations are formed as follows 

F Si+ 1 - F Si = sin(sj + i) - sin(sj), i = 0, ...,7V. (4.4) 

The fluxes, F Si , approximate solution values at the flux locations, and are computed by 
fully-upwind extrapolations from the primal nodes (except for the first interior dual node) as 


F s o — Uq, 

p _ (si—x 0 )U 1 + (xi — si)Uq 

Sl X±—Xo ’ 

p _ (Si—Xi- 2 )Ui-l-(.Si—Xi-l)Ui-2 

Si Xi-l—Xi-2 


i = 2, 


,N + 1 


(4.5) 


where (/,; is a discrete approximation to U(xi). These inviscid fluxes do not include solution 
derivatives and thus, according to (3.10), J _1 (Q) = 0(1) in all cases. All fluxes (4.5) are 
2 nd -order accurate and, thus, R ( Q ) = 0{h, 2 ). The analysis predicts 

Ed = 0(h 2 ) and E t = 0{h). (4.6) 


Figure 4.2 shows convergence rates obtained in grid-refinement and DS computations 
on irregular discretization grids involving random primal meshes and unbiased unperturbed 
dual meshes. The equivalent meshsize is taken as h = jj, where N is the number of primal 
intervals in grid-refinement computations. The convergence history of the L 0 c and L \ norms 
of truncation and discretization errors observed in the global computations confirms the sharp 
estimates obtained by the analysis (4.6) and in the DS tests. 

Table 4.1 summarizes discretization and truncation error convergence rates observed in 
computations with ID node-centered FVD schemes. Although not shown, similar conver- 
gence rates are observed for cell-centered FVD schemes. The results reveal two important 
trends valid for general non-degenerate inviscid equations: (1) grid irregularity strongly af- 
fects the truncation-error convergence, but has no effect on convergence of the discretiza- 
tion errors; (2) bias does not affect convergence rates. The analysis estimates are sharp for 
discretization errors on all grids and for truncation errors on all “randomized” grids. The 
convergence of truncation errors on uniform unperturbed grids is 2 nd order because of er- 
ror cancellations occurring on these regular grids. The DS-test convergence rates accurately 
predict the corresponding rates observed in global computations on all grids. Note that in 
multidimensional computations on regular grids, DS tests usually overestimate convergence 
of discretization errors. 
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(a) DS test (b) Grid refinement computations 


FIG. 4.2. Convergence of discretization and truncation errors in the DS test (a) and in grid-refinement com- 
putations (b) for the constant-coefficient convection equation. The tests are performed with random primal meshes 
and unbiased unperturbed dual meshes. 


Primal 

Dual 

DS test 

Grid-refinement computations 

Mesh 

Mesh 

Discr. Error 

Tmnc. Error 

Discr. Error 

Trunc. Error 

Uniform 

Unbiased, 

Unperturbed 

0{h“) 

0(h“) 

O(h') 

om 

Uniform 

Biased, 

Unperturbed 

om 

om 

0(h') 

0(h ‘) 

Uniform 

Unbiased 

Perturbed 

om 

o{h) 

O(h') 

O(h) 

Uniform 

Biased 

Perturbed 

om 

0{h) 

om 

o(h ) 

Random 

Unbiased, 

Unperturbed 

om 

0(h) 

om 

O(h) 

Random 

Biased, 

Unperturbed 

om 

0(h) 

om 

o(h ) 

Random 

Unbiased, 

Perturbed 

O(h') 

o{h) 

O(h') 

O(h) 

Random 

Biased, 

Perturbed 

om 

0(h) 

om 

o(h ) 


Table 4.1 

Convergence of discretization and truncation errors for node-centered FVD schemes for the constant- 
coefficient convection equation. 


4.3. Diffusion equation. The second set of one-dimensional tests illustrates the applica- 
tion of the DS analysis methodology to the diffusion equation. The non-degenerate constant- 
coefficient equation 


d xx U = f(x), U(0)=U o , U(1) = U U (4.7) 

is defined on the interval x £ [0, 1], with the exact solution U = sin(rr), / = — sin(s). FVD 
equations are formed as 

F Si+1 ~ F Si = cos(sj + i) - cos(si), i = l,...,N] U 0 = U 0 ; U N = U 1 . (4.8) 
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Fluxes approximating the solution derivative are defined as 


F Kt = 


Ui - Ui- 


i- 1 


i = 1, ...,7V, 


T>% %i—l 


(4.9) 


where U, is a discrete approximation to U(xi). Dirichlet boundary conditions are enforced. 
According to (3.10), J~ 1 (Q) = 0(h). 


Primal 

Dual 

; DS test 

| Grid-refinement computations | 

Mesh 

Mesh 

Discr. Error 

Trunc. Error 

Discr. Error 

Trunc. Error 

Uniform 

Unbiased. 

Unperturbed 

0(h *) 

om 

om 

om 

Uniform 

Biased, 

Unperturbed 

om 

0{h) 

o(h) 

O(h) 

Uniform 

Unbiased, 

Perturbed 

0(h?) 

0(1) 

om 

0(1) 

Uniform 

Biased. 

Perturbed 

om 

0(1) 

o(h) 

0(1) 

Random 

Unbiased. 

Unperturbed 

om 

0(h) 

om 

O(h) 

Random 

Biased, 

Unperturbed 

om 

0(1) 

o(h) 

0(1) 

Random 

Unbiased. 

Perturbed 

om 

0(1) 

om 

0(1) 

Random 

Biased. 

Perturbed 

om 

0(1) 

o(h) 

0(1) 


Table 4.2 


Convergence of discretization and truncation errors for node-centered FVD schemes for the diffusion equa- 
tion. The analytical estimates accurately predict DS-test convergence rates for all grid types except uniform- 
primal/unperturbed-dual type. The DS tests accurately predict convergence orders of grid-refinement truncation 
errors. The discretization-error convergence order is predicted well for the unbiased perturbed dual meshes, but 
overestimated in other cases. 


The placement of dual nodes (flux locations) significantly affects the error convergence. 
For unbiased, unperturbed dual meshes, the fluxes are approximated with the 2 -order ac- 
curacy, thus, providing R(Q) = 0(h 2 ); for such FVD schemes, the analysis predicts 

E d = 0(h 3 ) and E t = 0(h). (4.10) 

For either perturbed or biased dual meshes, R(Q) = 0(h), and the corresponding estimates 
are 


E d = 0(h 2 ) and E t = 0(1). (4.11) 

Table 4.2 summarizes the discretization and truncation error convergence rates observed 
for ID node-centered FVD schemes. The two main observations about grid-refined conver- 
gence rates of discretization errors are (1) local “random” grid irregularities do not affect the 
discretization-error convergence and (2) grid bias leads to one-order convergence deteriora- 
tion. 

The analytical estimates (4.10) and (4.11) accurately predict DS-test convergence rates 
for all grid types except uniform-primal/unperturbed-dual type; the convergence rates on 
these regular grids are faster because of local error cancellation. The DS tests accurately 
predict convergence orders of grid-refinement truncation errors. The discretization-error con- 
vergence order is predicted well for the unbiased perturbed dual meshes, but overestimated 
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in other cases. The reason for this overestimation is the error accumulation occurring in grid- 
refinement computations for all dual meshes except the unbiased perturbed type. Based on 
experience to date, in multidimensional computations on general irregular grids, the estimates 
obtained by the analysis and in DS tests are sharp as shown in following sections. 

5. Two-dimensional computations. 



FIG. 5. 1 . A typical mixed-element unstructured grid generated with random splitting and random perturbation 
of the underlying quadrilateral grid. 


5.1. Poisson equation. As a two-dimensional scalar example, we solve the Poisson 
equation. 


A U = f, (5.1) 

with Dirichlet boundary conditions, on a series of primal mixed-element grids composed 
of triangles and quadrangles. Each grid is formed from an underlying regular quadrilateral 
grid (Figure 5.1). In terms of a polar, (r, 9), coordinate system, the grid extent is defined 
as 6 £ [7t/6,7t/3] in the circumferential direction and r £ [1,2.2] in the radial direction. 
Grid irregularity is enforced by random splitting (or not splitting) quadrangles into triangles; 
approximately half of the quadrangles are split. In addition, the interior grid points are per- 
turbed from their original position by random shifts in the range [— -\/2/6, \/2/6] of the local 
mesh size in the radial direction. 


P 2 



FIG. 5.2. Illustration of gradient reconstruction for viscous terms on mixed grids with median-dual partition. 


The exact solution and forcing term are taken as U = [(sin(7rx)) 2 + (sin(7ry)) 2 ]/2, 
/ = — 27r 2 [l — (cos(7nr)) 2 — (cos(7 rj/)) 2 ]. In the FVD scheme, the solution is represented as a 
piecewise polynomial function, with polynomials defined at the primal cells; the conservation 
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law 


/ 


VU-ndr = 



(5.2) 


is enforced on node-centered control volumes constructed by the median-dual partition. 

With reference to Figure 5.2, the integral flux through the dual faces adjacent to the edge 
[Po, P4] is computed as 


/ VU-ndr = VU fl -n R + VU L n L , (5.3) 

ABC 

where n R and iil are directed areas of the corresponding dual faces. The gradient is recon- 
structed separately at each dual face as follows. 

For the triangular element contribution, the gradient is determined from a Green-Gauss 
evaluation at the primal-grid element, 

VU l =VUoi 4 . (5.4) 


The gradient overbar denotes a gradient evaluated by the Green-Gauss formula on the primal 
cell identified by the point subscripts. With fully-triangular elements, the formulation is 
equivalent to a Galerkin finite element scheme with a linear basis function [2, 3], For the 
quadrilateral element contribution, the gradient is evaluated as 


VU fi 


VU0234 + 


'V 4 
>4 


Po 

r ol 


— VU0234 • eo 4 


eo4, 


where r, is the coordinate vector of node P, and 


(5.5) 


eo4 


1~4 - fp 

| r 4 - r 0 | 


(5.6) 


is the unit vector aligned with the edge [Po, P4]. Note that for grids with dual faces perpen- 
dicular to the edges, the edge-gradient is the only contributor. This approach to the gradient 
reconstruction is used to decrease the scheme susceptibility to odd-even decoupling [11, 12], 
In all cases, the linear solution reconstruction leads to the I s * -order flux (gradient) recon- 
struction accuracy. 

The sequences of globally refined grids are generated with 2 n+3 + 1 points in both the 
radial and circumferential directions, where n = 0, 1, 2, 3, 4. The volume integral of the 
forcing term in (5.2) is evaluated as the value at the node multiplied by the dual volume. This 
integration is first order locally, but second-order convergence of this integration method over 
any 0(1) domain is observed on general mixed-element grids. The sequences of DS grids are 
generated from a grid with 17 points in both the nominal radial and circumferential directions 
and downscaled by the scaled-grid method about the center of the domain by a factor 2~ s , 
where s = 0, 2, 4, 6, 8. The grid topology remains unchanged. 

The asymptotic order of computed according to (3.10), is m j = 1. Overall, 

R(Q) = 0(h) and the analysis predicts 

E d = 0(h 2 ) and E t = 0(1). (5.7) 


The L 1 norms of truncation and discretization errors are shown in Figure 5.3 versus an 
equivalent mesh size parameter, taken as the Li norm of the square root of the dual volume. 


15 




- 1st order slope 

2nd order slope 

© Discretizaton Error LI -norm 

B Truncation Error LI -norm 



i I i i i I i I i i i I, 

0.02 0.04 0.06 0.08 

Equivalent Mesh Size 


(a) DS test 


(b) Grid refinement computations 


FIG. 5.3. Convergence of the discretization and truncation errors for the Poisson equation solved on mixed- 
element unstructured grids. 


Although not shown, error convergence rates in the norm are the same as the L \ -norm 
rates. The analysis (5.7) and the DS tests predict grid-refinement error convergence rate 
precisely; the truncation errors remain 0(1) and the discretization errors converge with 2 nd - 
order. The reason for the 0(1) convergence of truncation errors is grid irregularity stemming 
from the perturbation to the grid points and the usage of mixed grids. Both references [19] 
and [22] interpret 0(1) convergence of truncation errors on irregular grids as indication that 
the corresponding discrete solutions do not approximate the continuous ones; this example 
clearly shows that this is not the case. Although not shown, with regular meshes composed of 
either triangles alone or quadrangles alone, both the truncation errors and the discretization 
errors converge with 2 nd order. 

5.2. Incompressible Euler equations. In this section we consider incompressible in- 
viscid flow equations. The source function S is assumed to be zero. Inviscid fluxes for 
conservation of mass and momentum are defined as 



flu 


flv 

F = f i + gj = 

u 2 +p 

i + 

uv 


uv 


v 2 +p 


where the vector of unknowns, q = [m, v,p], includes the Cartesian velocities and the pres- 
sure; fl is an artificial compressibility parameter introduced as in [2] and taken as fl = 1 
here. 

The median-dual partition is applied. At each dual control volume, a polynomial solution 
approximation is constructed. The approximation is required to coincide with the discrete 
solution value at the central node Pq\ the polynomial coefficients are determined through 
a least-square procedure involving neighboring nodes, i.e., nodes linked by an edge to the 
central node /),. If the set of neighboring nodes is insufficient to determine, uniquely and sta- 
bly, the polynomial coefficients, the set may be expanded to involve neighbors of neighbors. 
Three node-centered FVD schemes are considered: an edge-reconstruction scheme and two 
face-reconstruction schemes. 
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5.2.1. Equivalent linear operator analysis. The Jacobian of <f (F • n) dT is a 3 x 3 

r 


matrix 


Tly 0 \ 

un y n x dT (5.9) 

un x + 2vn v n y J 

operating on the vector 5q = [du, 5v, 5p] T . After equivalent substitutions, the equivalent 
linear operator, E(Q), becomes 

( 0 ( 1 ) 0 ( 1 ) 0 \ 

E(Q)=\ uO{l) + uO(l) uO{ 1) 0(1) . (5.10) 

\ vO( 1) uO(l) + vO(l) 0(1) / 

If u 2 + v 2 > e > 0 everywhere in the computational domain, i.e., the system is non- 
degenerate, the inverse of E(Q) is a full matrix 

/ 0 ( 1 ) 0 ( 1 ) 0 ( 1 ) \ 

J~ 1 {Q) ~ E~ 1 (Q) = 0(1) 0(1) 0(1) . (5.11) 

V 0(1) 0(1) 0(1) J 

Thus, the estimate of the discretization error convergence is fully determined by convergence 
of contributors to the residual R{Q) (Section 3.4.1). In a general situation, the asymptotic 
order of J~ 1 (Q) can be different for different equations, but for non-degenerate (away from 
stagnation) inviscid flow equations, the asymptotic order is mj = 0 for all three equations. 
The effect of degeneration (flow stagnation) on convergence of the discretization errors is con- 
sidered subsequently in Section 5.3. In the course of Section 5.2, we assume non-degenerate 
equations implying non-vanishing velocity components of the (manufactured) solution. 




FIG. 5.4. Illustration for edge-reconstruction flux integration scheme in the interior. 


5.2.2. Edge-reconstruction scheme. This edge-reconstruction FVD scheme is widely 
used in unstructured-grid computations [2, 3, 11, 12, 19]. The term edge-reconstruction 
emphasizes that the quadrature scheme used for computing the integrals over the dual-cell 
boundaries employs fluxes reconstructed at the midpoints of the edges connecting the grid 
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nodes. By construction, the dual-cell boundaries are piecewise linear and, therefore, in the 
interior and next to straight physical boundaries, the outward normal n approximation is ex- 
act. 

The numerical upwind flux at the edge midpoint is approximated according to the Roe 
scheme [2, 21], 

(F h -n)=l [(F 0 • n) + (F, • n)] - ^ |A(Q)| (Q x - Q 0 ), (5.12) 

where, Qo and Qi are the “left” and “right” solution reconstructions at the edge midpoint 
derived from the linear approximations at the control volumes centered at Pq and Pi, respec- 
tively (see Figure 5.4); Fo and Fi are the corresponding “left” and “right” numerical fluxes; 
n is the unit vector aligned with the combined-directed-area vector n = iil + rip , where 
riL and iir are outward normal (directed-area) vectors of the left and right segments, respec- 
tively, with their lengths equal to the corresponding segment lengths; | A(Q)| is the Roe’s 
approximate Riemann solver matrix computed for Q = -) ( Q i + Qo). With a linear solution 
approximation at each control volume, this is a 2 -order accurate flux reconstruction. 

The flux integration over the two segments of the control-volume boundary linked at 
the edge midpoint is approximated by multiplying F h computed at the edge midpoint with 
the length of the combined-directed-area vector, |n|. This computationally efficient scheme 
provides exact integration for the conservation laws with linear fluxes on triangular grids, 
where it corresponds to a linear finite-element Galerkin scheme [2, 3, 4], According to the 
definition introduced in the sufficient flux-integration accuracy (condition 3 in Section 3.4.1), 
the scheme possesses global 2 -order accuracy. The scheme provides only 1 sf -order local 
integration accuracy. On general (irregular) quadrilateral and mixed-element grids, the global 
flux integration accuracy deteriorates to 1 st order; examples are shown subsequently in Sec- 
tion 5.2.4; additional examples can be found in [8], 

The forcing term integration over the control volume is approximated as the node value 
multiplied by the volume |fi|. This integration method is first-order accurate locally, but 
the integration over any 0(1) domain is second order on general globally refined grids with 
a median-dual partition. Boundary conditions are enforced weakly through the boundary 
fluxes; implementation details are given in Appendix A. The overall residual convergence 
is the same for all equations and is estimated as R(Q) = 0(h 2 ) on triangular grids and 
R(Q) = 0(h) on general quadrilateral and mixed-element grids, leading to estimates 


Ed, = 0(h 2 ) and E t = 0(h ) 

(5.13) 

for irregular triangular grids and 


Ed = 0{h ) and E t = 0(1) 

(5.14) 


for general irregular grids. 

5.2.3. Face-reconstruction schemes. In this section, we describe two face-reconstruc- 
tion FVD schemes that employ the median-dual partition, and provide 2 nd and 3 rd order ac- 
curacies on general irregular grids. The term face-reconstruction refers to the Gauss-Legendre 
method for flux integration; the method uses fluxes reconstructed at Gaussian points defined 
at each segment (face) of the dual-cell boundaries. Similar to the edge-reconstruction scheme, 
the linear and quadratic polynomials are defined at dual cells and coincide with the solutions 
at the grid nodes. The polynomial coefficients are defined in a least-square procedure involv- 
ing solutions at the neighboring nodes. In the interior, piecewise-straight dual boundaries 
imply an exact n. The schemes described in this section are quite similar to the scheme dis- 
cussed by Delanaye and Liu [6] for cell-centered discretizations. As noted therein and shown 
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in [8], the operation counts significantly favor the cell-centered approach for flux integration 
in three dimensions. 

The distinguishing feature of these face-reconstruction schemes is the flux integration 
procedures ensuring the designed local integration accuracy. Extension of 2 -order accu- 
racy to mixed grids can be achieved with linear polynomials, a modified straight-segment 
tangency boundary approximation (discussed subsequently in Section 5.2.4), and a flux inte- 
gration scheme providing local 2 -order accuracy. Accuracy of 3 rd order requires quadratic 
polynomials for flux reconstruction, quadratic fit to the curved tangency boundary, and an 
flux integration scheme with local 3 -order accuracy. 

In the implemented 2D version, a polynomial flux is defined at each segment of the dual 
control-volume boundary and used in (5.12). The “left” and “right” solutions, Qo and Qi, 
are represented by the basic polynomials defined at the adjacent control volumes; numerical 
fluxes, Fo and Fi, are analytically computed as products of the corresponding basic poly- 
nomials. The dissipation matrix, | A(Q)|, is computed for Q, defined as the average of Qo 
and Qi evaluated at the segment midpoint. The polynomial flux at this segment is defined 
according to (5.12), where all terms, beside | A(Q)|, are polynomials. The scheme with lin- 
ear polynomials provides 2 -order accurate flux reconstruction; fluxes reconstructed with 
quadratic polynomials are 3 rd -order accurate. 

The inner product of the polynomial flux vector and the outward unit normal vector is 
integrated over each segment of the dual-cell boundary with a numerical Gauss-Legendre 
quadrature formula employing 1 point for the linear face-reconstruction scheme, and 2 points 
for the quadratic face-reconstruction scheme. As a convenient debugging tool, one can add 
one Gaussian point to have total of 2 and 3 points, respectively, per integration segment. 
With these improved accuracy integrations, the FVD method should provide zero residuals 
for linear and quadratic manufactured solutions. 

The forcing term integration method for the linear face-reconstruction scheme is the 
same as for the edge-reconstruction scheme, providing 2 nd -order accuracy. For quadratic 
face-reconstruction scheme, a linear (or piecewise linear) approximation to the forcing term 
is constructed and integrated at each control volume providing a 3 rd -order accurate approx- 
imation to Tjn- ff f dCl. Description of boundary condition implementation for the face- 
ts 

reconstruction FVD schemes is provided in Appendix A. Overall, R{Q) = 0{h 2 ) for the 
linear scheme and R(Q) = 0{h 3 ) for the quadratic scheme; residual convergence does not 
deteriorate for general irregular grids. The predictions for the error convergence rates are 

E d = 0(h 2 ) and E t = 0{h ) (5.15) 

for the linear FVD scheme and 

E d = 0(h 3 ) and E t = 0(h 2 ) (5.16) 

for the quadratic FVD scheme. 

5.2.4. Numerical tests for non-degenerate flows. Numerical tests presented in this 
section are performed for 2D inviscid incompressible flows around a cylinder of unit radius 
centered in the origin. The flow is described by the conservative equations (1.1) with zero 
source and forcing terms and fluxes defined in (5.8). The analytical solution for this problem 
is known 


U = Uoo + 


2 sin 2 6-1 


V = Voo + 2^ 

p _ p _ u 2 +v 2 

r — r oo 2 
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+ dylp, 

- d x !p, 


(5.17) 



where (r, 9) are the polar coordinates r 2 = x 2 + y 2 , tan 9 = y/x, and -0 = — Cln(r) is the 
stream function with C being a constant characterizing the flow circulation. In the course of 
this section, the free stream at the infinity is characterized by Uoo = 1, Voo = 0, = 1.5, 

and the zero circulation (C = 0) is assumed. 





(a) Irregular triangular grid (b) Irregular quadrilateral grid (c) Irregular mixed grid 


FIG. 5.5. Typical triangular and mixed-element unstructured grids defined on a computational domain shifted 
away from the surface of the cylinder. 


Interior inflow/outflow domains. The first set of tests is performed on a computational 
domain shifted away from the surface of the cylinder: 1.5 < r < 4, 2ir / 3 < 9 < 47t/3. 
Seven formulations are studied: the edge-reconstruction FVD scheme on irregular triangular, 
quadrilateral, and mixed-element grids; and face-reconstruction FVD schemes on irregular 
quadrilateral and mixed-element grids. Examples of irregular grids derived from an underly- 
ing regular grid are shown in Figure 5.5. For triangular and mixed-element grids, irregularity 
is introduced through random splitting (or not splitting) of structured quadrilateral cells. Each 
cell has equal probabilities to introduce either of the two diagonal choices or, for mixed- 
element grids, no diagonals. For irregular quadrilateral grids, interior nodes are perturbed 
from their original position by random shifts in the range [ — \/2/6, \/2/6] of the local mesh 
size in both the radial and circumferential directions. 

For each formulation, grid refinement and DS tests are performed. In grid refinements, 
the underlined regular grid is refined by doubling the number of intervals in the radial and 
angular directions. In the DS test, the coarsest 9x9 grid is scaled down around the point 
r = 2.75, 9 = 7r by multiplying all angular and radial deviations from this point by a factor 
of 0.5. Randomization is introduced independently on each scale. The inflow boundary con- 
ditions are enforced at the boundary corresponding to the external radius; outflow conditions 
are enforced at all other boundaries. Table 5.1 summarizes the convergence of discretization 
and truncation errors observed in these tests. The results confirm the analytical predictions 
(5.13) — (5.16) and the capabilities of the DS test to provide sharp estimates for error con- 
vergence in grid-refinement computations. 

The convergence rates observed for the edge-reconstruction scheme on irregular triangu- 
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Formulation 

Analysis / DS test 

Grid-refinement 

computations 

Discr. Error 

Trunc. Error 

Discr. Error 

Trunc. Error 

Edge-reconstruction, 
irregular triangular grid 

0(h 2 ) 

O(h) 

0{h 2 ) 

0(h) 

Edge-reconstruction, 
irregular quadrilateral grid 

o(h ) 

0(1) 

0(h) 

0(1) 

Edge-reconstruction, 
mixed-element grid 

0(h) 

0(1) 

0(h) 

0(1) 

Linear 

face-reconstruction, 
irregular quadrilateral grid 

0(h 2 ) 

O(h) 

o(h 2 ) 

0(h) 

Linear 

face-reconstruction, 
mixed-element grid 

o(h 2 ) 

0(h) 

o(h 2 ) 

O(h) 

Quadratic 

face-reconstruction, 
irregular quadrilateral grid 

o(h 3 ) 

o(h 2 ) 

o(h 3 ) 

o(h 2 ) 

Quadratic 

face-reconstruction, 
mixed-element grid 

o(h 3 ) 

o(h 2 ) 

o(h 3 ) 

o(h 2 ) 


Table 5.1 


Convergence of discretization and truncation errors for various irregular grid formulation of the 2D inviscid 
incompressible equations on an inflow/outflow computational domain. The error convergence orders predicted by 
the analysis are observed for all variables and equations in all norms. 


lar and irregular quadrilateral grids are consistent with the results reported in [1], Although 
not shown, we have implemented a central version of the edge-reconstruction scheme, where 
fluxes at the edge midpoints are defined as averages of the solutions in the neighboring nodes 
The observed convergence rates were identical with the rates shown in Table 5.1. The re- 
sults contradict to [19], where zeroth-order convergence of discretization errors on irregular 
quadrilateral grids was reported for a central scheme for a constant-coefficient convection 
equation. An in-depth investigation of this discrepancy confirming correctness of our results 
has been reported in [8], 

Tangency boundary. For the edge-reconstruction FVD scheme on triangular grids, local 
accuracy deterioration occurs if a curved tangency boundary is approximated by straight seg- 
ments linking primal-mesh nodes located at the physical boundary. Sketch (a) of Figure 5.6 il- 
lustrates this approximation: the straight segments [P2 , Po] and [Po , P4] approximate a curved 
boundary, the points Pb and Pc are the segments’ midpoints, and the arrow indicates the lo- 
cal flow velocity. For the edge-reconstruction scheme, discrete tangency is enforced over the 
straight segments [Po, Pb] and [Po, Pc]- The exact continuous solution satisfies the tangency 
condition at the actual curved boundary, not at the straight boundary segments. If evaluated 
with the exact solution, one boundary segment would appear as inflow, while the other would 
appear as outflow. 

To assess the error introduced by implementing zero velocities over the straight boundary 
segments, consider the segment [Pb , Po] that contributes to the control volume centered at 
Po. An overspecified version of the boundary conditions, in which the exact velocities are 
reconstructed at P = 5/6Po + I/6P2, would complement the interior edge-reconstruction 
formulation to provide a 2"' / -order accurate solution [3], Thus, the difference between this 
overspecified boundary flux and the weak tangency boundary condition can be considered 
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(a) Basic median-dual approximation 
for the edge-reconstruction scheme 


(b) Modified approximation for the 
linear face-reconstruction scheme 


FIG. 5.6. Straight-segment approximations to curved tangency boundary. 


as a velocity reconstruction error. In particular, for the mass conservation law, the leading 
contribution to R[Q) from the segment [Pq, Pb\ can be estimated as 

R[Pb,Po\ = V • n = 0(h), (5.18) 

where 

V = |V 0 + ±V 2 = § (V 0 + V 2 ) + 0(h), 
n = \ (n 0 + n 2 ) + 0{h 2 ). 

and 


(V 0 + V 2 ) • (n 0 + n 2 ) = 0(h 2 ), 

Here, V is the velocity vector reconstructed at P and ii is the outward unit normal correspond- 
ing to [Pq, Pb\. Vo, no, V 2 , n 2 , are the exact velocities and unit normals at the grid nodes 
Po and P 2 , respectively, satisfying (Vo • no) = (V 2 • n 2 ) = 0, and h is the characteristic 
meshsize. The l st -order accuracy in R(Q) leads to local 1 st -ordcr accurate discretization er- 
rors. The reason for the residual contribution, R[p B t p 0 i, to be 0(h) is the l st -order difference 
between V and (Vo + V 2 )/2. Similar 0(h) contributions appear in momentum-conservation 
equations. 

To illustrate the effect of a straight-segment approximation to a curved tangency bound- 
ary, a sequence of irregular triangular grids is generated at the top of the cylinder (1 < r < 
2.2, 7t/3 < 9 < 27t/ 3) and used in computations with the edge-reconstruction FVD scheme; 
a grid example is shown in Figure 5.7. Figure 5.8 illustrates convergence of the L \ norm 
of truncation and discretization errors in grid-refinement and DS tests performed with edge- 
reconstruction FVD scheme. Two DS tests are performed, each with the focal point at the top 
surface of the cylinder. The first DS test uses overspecification at all boundaries except the in- 
terior tangency nodes. The second DS test replaces the overspecification along one boundary 
with the physical inflow boundary condition (see sketches in Figure 5.8). 

For the interior-tangency nodes, there are two tangency segments; the errors at these seg- 
ments contribute to R(Q) with opposite signs and, at least partially, compensate each other. 
Because of this compensation, the accuracy deterioration does not affect interior-tangency 
nodes on the grids with (nearly) uniform boundary node distributions. The 2 -order conver- 
gence of discretization errors and the l st -order convergence of truncation errors demonstrated 
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1.5b 


FIG. 5.7. Irregular triangular grid around the top of the cylinder. 



(a) DS test: interior tangency bound- 
ary conditions 



(b) DS test: inflow/tangency boundary 
conditions 


10' w 



Equivalent Mesh Size 



1st order slope 
2nd order slope 
Truncation Error : Continuity 
Truncation Error : x- momentum 
Truncation Error : y- momentum 
Discretization Error : p 
Discretization Error :u 
Discretization Enor : v 


(c) Grid refinement 

FIG. 5.8. Convergence of the L± norm of truncation and discretization errors observed in DS and grid- 
refinement tests for the edge-reconstruction FVD scheme. The tests are performed on irregular triangular grids 
surrounding the top tangency boundary of the unit cylinder. The open square symbols in the sketches denote over- 
specification for the DS tests. 
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in the interior-tangency DS test confirm this conclusion. However, in the corners or/and at the 
interfaces between clusters of boundary nodes with different topology/metrics, the compen- 
sation does not occur. The accuracy deterioration is clearly observed in the DS test performed 
with the inflow/tangency boundary conditions. The convergence of the L i norms of the errors 
in the grid-refinement test is the same as convergence shown in the interior-tangency DS test 
and is not affected by local accuracy deterioration in the corners; although not shown, the 
grid-refinement convergence of the L 0 0 norm of errors is slower and similar to the conver- 
gence demonstrated in the DS test with the inflow/tangency boundary conditions. 




(a) DS test: basic approximation 


(b) DS test: modified approximation 



1st order slope 
2nd order slope 
Truncation Error : Continuity 
Truncation Error : x- momentum 
Truncation Error : y- momentum 
Discretization Error : p 
Discretization Error : u 
Discretization Error : v 


FIG. 5.9. Convergence of the L\ norm of truncation and discretization errors observed in DS tests for the 
linear face-reconstruction FVD scheme with two types of straight-segment approximation to the curved tangency 
boundary. 


The linear face-reconstruction scheme possesses the flexibility to recover the 2 -order 
accuracy with a straight-segment tangency boundary approximation. Two required modifica- 
tions are illustrated in the sketch (b) of Figure 5.6: (1) the points Pb and Pc are moved to 
the boundary, and (2) the fluxes are reconstructed at the midpoints of the straight boundary 
segments, e.g., at P = 1/2Pb + 1/2Pq. Figure 5.9 shows convergence of the L\ norm of trun- 
cation and discretization errors observed in DS tests for the linear face-reconstruction FVD 
scheme with two types of straight-segment approximation to the curved tangency bound- 
ary: the basic median-dual approximation and the modified approximation. The tests are 
performed in the setting similar to the inflow/tangency boundary conditions used with the 
edge-reconstruction scheme. The results confirm that, with modifications (1) and (2), the 
2"' / -order convergence can be achieved with linear polynomials and a straight-segment ap- 
proximation to the curved tangency boundary. Slower convergence observed with the basic 
median-dual cells indicates that both modifications are essential. Although not shown, with 
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the quadratic face-reconstruction scheme and a quadratic fit to the curved tangency boundary, 
discretization errors converge with the 3 rd order. 

5.3. Discretization error deterioration at stagnation. Local discretization accuracy 
deterioration may occur in the vicinity of stagnation or near other singularities in degener- 
ated equations. To provide an insight into this phenomenon, we consider a one-dimensional 
symmetry-line nonlinear convection equation 


uu x = f 


(5.19) 


linearized around the leading-edge stagnation solution for the flat-plate geometry ( u = —x) 


-xd h x b u — S w = g, 6 U (- 1) = (5i nflow , (5.20) 

where d x is a discrete derivative approximation; the same equation describes the tangency 
flow in the vicinity of the trailing-edge stagnation. The solution of (5.20) represents the 
discretization error, S u . The right-hand side g is 0(h p ) small for p th -order discretizations and 
represents possible residual perturbations; the inflow condition, S u = 0 defined at a: = — 1, 
may also contain an <5i n fl ow = 0(h p ) error. The general solution of the initial-value problem 
(4.2) on the interval x £ [0, 1] is given by 


A" — A" f'r'l 

u ^particular v**' / 


■^particular (-1)- ^inflow 


(5.21) 


where <5p articular (a:) is a particular solution of the equation (5.20). The function ^p art i cu i a r( ;c ) 
and the constant (5; n fl ow are both 0(h p ) small. However, the discretization error S u grows as 
0(1 /x) in the vicinity of stagnation ( x = 0) and locally becomes 0(/i p_1 ). 

Analysis of the equivalent linear operator performed for the inviscid equations is also 
capable of detecting the discretization-error deterioration. In the vicinity of stagnation, both 
velocity components become Oft) small and the equivalent operator becomes 


/ 0 ( 1 ) 0 ( 1 ) 0 \ 

E{Q) = 0(h) 0(h) 0(1) , 

V 0(h) O(h) 0(1) ) 

cf. (5.10). Thus, 


/ 0(1) Oft" 1 ) Oft- 1 ) \ 

J~\Q) ~ E~\Q) = 0(1) Oft" 1 ) Oft" 1 ) 

v Oft) 0(1) 0(1) ) 


(5.22) 


(5.23) 


For the velocity discretization errors, the asymptotic order of J~ X (Q) becomes mj = 
— 1, but for the discretization error in pressure, the asymptotic order remains the same as 
in the interior, rnj = 0, implying different orders of convergence for different variables at 
stagnation. To retain uniformly the same discretization-error convergence for all variables, 
one has to approximate the momentum conservation equations at stagnation with higher order 
than in the interior. 

Discretization accuracy deterioration in the vicinity of stagnation is a universal phe- 
nomenon and expected for any discretization of the inviscid flow equations. In simulations. 
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the stagnation accuracy deterioration is rarely noticed and seldom correctly attributed because 
its effect is mostly visible far downstream of stagnation. At stagnation, the solution and all 
its derivatives are very small; therefore, in spite of increased relative discretization error, the 
absolute value of the local error, which is proportional to some derivative of the solution, does 
not exceed a typical error away from stagnation. Also, the effect disappears in the presence 
of the diffusion terms that prevent equation degeneration. 

The discretization error generated at stagnation is amplified by the accelerated flow con- 
vecting the error downstream. For leading-edge stagnation, this error is convected along the 
tangency boundary, affects very few points, and is almost not observable in integral norms, 
such as L i and L 2 norms, but can be clearly detected in the L 0 c norm. At the trailing-edge 
stagnation, the generated discretization error is convected into the interior affecting a larger 
area. This may explain the difficulties observed in solution of the trailing edge flow problems, 
even for relatively benign geometries such as flat-plate, cylinder, or parabola. 

In order to see accuracy order reduction more clearly, the DS test at stagnation needs 
to be adjusted to evaluate the relative discretization error, defined as \Ed\ / \Q\, where, for 
computing the relative discretization error in the velocities, \Q\ vel = max (Vu 2 + v 2 ) and, 
for the relative discretization error in the pressure, \Q\ p = max (p); maxima are taken over 
the DS-test domain. In fact, evaluation of discretization error is always required in the form 
of relative error, independent on the solution magnitude. Away from stagnation, the absolute 
and relative errors converge with the same orders because \Q\ vet = 0(1) and \Q\ = 0(1); 
at stagnation, \Q\ vel = 0(h) and \Q\ p = 0(1). 


Formulation 

| Error convergence at trailing-edge stagnation, O (h a ) j 

DS test 

] Grid-refinement computations j 

Relative 
Discr. Error 

Trunc. 

Error 

| Discr. Error 

Trunc. 

Error 

Li norm 

Loo norm 

Edge-reconstruction, 
irregular triangular grid 

u, v : a = 1 
p : a — 2 

a = 1 

1 < a < 2 

a = 1 

a = 1 

Edge-reconstruction, 
mixed-element grid 

u, v : a = 0 
p : a = 1 

a = 0 

0 < a < 1 

a = 0 

a = 0 

Linear 

face-reconstruction, 
mixed-element grid 

u, v : a = 1 
p : a = 2 

a = 1 

1 < a < 2 

u,v : a = 1 
p : 1 < a < 2 


Quadratic 

face-reconstruction, 
mixed-element grid 

u, v : a = 2 
p : a — 3 

a — 2 

2 < a < 3 

a = 2 



Table 5.2 


Convergence of discretization and truncation errors for various irregular grid formulations of the 2D inviscid 
incompressible equations at the aft of the unit cylinder. A single entry of a in the table indicates it refers to all the 
variables (u, v, p). One order reduction in convergence of discretization errors for velocity components is clearly 
observed in DS tests and in the Loo norm of the discretization errors obtained in grid refinement (cf Table 5.1); 
slowdown of L i -norm convergence is also observed. Truncation errors converge with the same order as in the 
interior. 


Table 5.2 summarizes the error convergence orders observed in computations at the aft 
of the unit cylinder. Each global grid is formed from an underlying regular quadrilateral grid 
spanning 120 degrees in 9 (centered about the most downstream point on the cylinder surface) 
with extent 1 < r < 2.2 in the radial direction; the grids are generated with 2" +3 + 1 points 
in both the radial and circumferential directions, where n — 0, 1, 2, 3, 4. On each grid, the 
regular quadrangles are split randomly into a fully-triangular grid or a mixed-element grid, 
as discussed previously. The DS test is performed with independent grid generation and the 
focal point placed at the most downstream point on the cylinder. 
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The analytical estimates and the DS test accurately predict convergence order of trun- 
cation errors and one-order degradation of the L 0 c norm of discretization errors in the grid- 
refinement computations. The integral L \ norm of the discretization error is less sensitive to 
the degradation; the error generated at stagnation is propagated along the streamline coming 
from the stagnation, and the interior areas far away from the streamline are little affected. 
Note also that a clear distinction between convergence of relative discretization errors in the 
velocities and pressure predicted by (5.23) is observed only in the DS tests; in the global 
grid refinement computations, the nonlinear interactions slow the convergence of discretiza- 
tion errors in pressure, although, the pressure errors remain much smaller than the velocity 
discretization errors on the same grids. 

6. Discussion. A new computational analysis tool, downscaling (DS) test, has been in- 
troduced and applied for studying the convergence rates of truncation and discretization errors 
of finite-volume discretization (FVD) schemes on general irregular grids. The study shows 
that the design-order convergence of discretization errors can be achieved even when trunca- 
tion errors exhibit a lower-order convergence or, in some cases, do not converge at all. The 
DS test is a general, efficient, accurate and practical tool, enabling straightforward exten- 
sion of verification and validation to general unstructured grid formulations. It also allows 
separate analysis of the interior, boundaries, and singularities that could be useful even in 
structured-grid settings. 

There are several new findings arising from the use of the DS test analysis. It was 
shown that the discretization accuracy of a common node-centered edge-reconstruction FVD 
scheme, known to be 2"' , -order accurate for inviscid equations on triangular grids, degener- 
ates to I s * order for mixed grids. Alternative node-centered face-reconstruction schemes have 
been presented and demonstrated to provide 2 nd and 3 rd order accuracies on mixed grids. The 
local accuracy deterioration at intersections of tangency and inflow/outflow boundaries has 
been demonstrated using DS tests tailored to examine the local behavior of the boundary con- 
ditions. The discretization-error order reduction within inviscid stagnation regions has been 
demonstrated. The accuracy deterioration is local, affecting mainly the velocity components, 
but applies to any order scheme. The result is somewhat surprising because the solution is 
so simple but analysis of the Jacobian operator along the stagnation streamline has provided 
insight into the phenomena. 

In 2-D inviscid computations, the cost of computing residuals of the face-reconstruction 
discretizations is about twice as large as the cost of edge-reconstruction residuals because the 
face-reconstruction discretizations require a solution of the approximate Riemann problem at 
each dual-cell boundary segment (face); the edge-reconstruction discretization requires one 
Riemann solution per two connecting segments. The cost increase in 3D computations for 
general tetrahedral grids is much larger because multiple dual control-volume faces are ad- 
jacent to each primal-mesh edge. A dramatic cost reduction for face-reconstruction schemes 
is achieved by evaluating the dissipation at the midpoints of primal-mesh edges, in the same 
way as in the edge-reconstruction discretization, once for all adjacent dual segments/faces; 
unsplit flux contributions are still evaluated at the dual segments/faces. Also, as shown in [8], 
the cost of face-reconstruction FVD schemes on cell-centered grids is much lower. 

On typical 2-D unstructured grids, where most nodes have at least five neighbors, a 
typical computational stencil for the quadratic face-reconstruction discretization has the same 
size as a stencil for the linear face-reconstruction scheme. The memory requirements for the 
quadratic discretization are about two times higher than for linear discretizations because a 
larger number of coefficients is stored at each control volume. 
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Appendix A. Boundary conditions for the incompressible Euler FVD schemes. 

For the edge-reconstruction FVD scheme near a curved boundary, the computational 
boundary is constructed as a piecewise straight line connecting the grid nodes located at the 
physical boundary. The control volume around a boundary node, I\ h (see Figure A.l) is 
closed with the boundary segments [Pb,Pq\ and [Pq,Pc\- The straight-line approximation 
[Po, P 2 ] provides a 2 -order accuracy to the curved boundary segments connecting nodes 
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FIG. A. 1. Illustration for edge-reconstruction flux integration scheme near the boundary. 


Po and P‘2- However, the approximation [Po, Pb } provides only a l st -order accuracy for the 
Po , Pf part of the curved boundary segment. The flux reconstruction and integration over the 
segments [Pa, Pe\ and [Pe, Pd] are the same as for the interior edges. The integral fluxes 
through segments [ Pa,Pb\ and [Pd, Pc] are approximated by the numerical fluxes ( 5 . 12 ) 
evaluated at points Pb and Pc using the directed areas ni and 114, respectively. The linear 
solution approximation at control volumes ensures 2 -order flux reconstruction accuracy. 

Boundary conditions are enforced weakly through the boundary fluxes. At the inflow, 
where velocity components are specified through prescribed boundary conditions, the bound- 
ary flux uses the specified velocities, the pressure is reconstructed from the discrete solution 
approximation. The outflow conditions are implemented similarly: the pressure contribution 
to the boundary flux is specified from prescribed conditions; the velocity components are re- 
constructed from the discrete solution. At the tangency boundary, the mass fluxes and velocity 
contributions to the momentum fluxes are explicitly set to zero; the pressure is reconstructed 
from the discrete solution. 

On triangular grids, the integral flux through the boundary segment [Pb , Ho] is computed 
as (F h [p B p 0 ] • 112), where 112 is the directed area of [Pb, Po], and F h [p s p 0 ] is evaluated at 
P = 5/6Po + I/ 6 P 2 . The solution components prescribed in the boundary conditions, e.g., 
both velocities at inflow, or the pressure at outflow, are specified at P from the known exact 
solution; other components are interpolated to P from the endpoints of the segment [P2, Po]. 
The integral flux through [Po, Pc] is computed analogously. In spite of providing only 1 st - 
order local accuracy for integrated fluxes through [Pb, Po] and [Po, Pc], the coefficients 5/6 
and 1/6 lead to cancellation of l st -order errors, providing zero residuals for conservation 
laws with linear fluxes and, thus, supporting the global 2 -order flux integration accuracy 

[3]. 

On mixed-element and irregular-quadrilateral grids, the 2 -order flux integration accu- 
racy is not recovered. For triangular cells adjacent to the boundary, the ( 5 / 6 , 1 / 6 ) rule is 
still enforced; on the boundary segments of quadrilateral cells, the fluxes are evaluated at the 
segment midpoint, P = 3/4Po + The overall boundary approximation accuracy for 

R(Q) is consistent with the accuracy of the interior FVD scheme, 2 nd order for triangular 
grids and 1 st order for other irregular grids. 

Reliance on error cancellation in providing design-order discretization accuracy may lead 
to some counterintuitive phenomena. In particular, improving the accuracy of boundary con- 
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ditions in a way that is not compatible with the flux computations through interior segments 
can upset the error cancellation balance and, in fact, worsen approximation accuracy of R(Q), 
at least, locally. In some cases, this local imbalance is compensated by contributions from 
other boundary segments. The imbalance is especially prominent at the corners where the 
boundary segments represent different types (e.g., tangency and inflow) of boundaries. 

For control volumes at general curved tangency boundaries, the tangency condition is 
enforced over the straight segments, rather than over the physical curved boundary. The er- 
ror introduced by this approach can be considered as a flux reconstruction error introduced 
by enforcing zero velocity over non-tangency boundary segments. As estimated analytically 
and confirmed numerically in Section 5.2.4, the error is 0(h) for each segment of the tan- 
gency boundary. Reference [16] investigates an alternative approach, in which other types 
of boundary conditions are enforced over the straight boundary segments approximating a 
curved tangency boundary. 

The boundary conditions for the face-reconstruction discretizations are also enforced 
weakly, through boundary fluxes. At inflow and outflow, the computational domain has a 
piecewise straight boundary; the solution components are either specified from the exact 
solution (the velocity is specified at inflow and the pressure is specified at outflow) or repre- 
sented by the polynomial approximation defined at the adjacent control volume. At tangency, 
the 2-order accuracy can still be achieved with enforcing tangency over piecewise linear 
segments. (See discussion in Section 5.2.4.) For the 3 rd -order accuracy, the physical bound- 
ary should be approximated quadratic ally. In a general case, where the analytical shape of the 
boundary is unknown, the boundary should be represented as a piecewise polynomial curve 
providing the required accuracy for the boundary shape. 
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